 use "${tmp}/Regfile_nafZE", replace

	
	xtset gp year
	
	
	destring ZE2010, replace
	
	// define sample 
	egen mlabexid = group(mainlabex)
	qui su mlabexid
	replace mlabexid = r(max) + ZE2010 if mi(mlabexid)
	

	reghdfe ln_eng_sbrut E_post11  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(ZE2010#year  gp)  cl(mlabexid)	
	bys gp: egen insampeng = max(e(sample))
	keep if insampeng
	 
	
	// Figure B2a	
	reghdfe ln_eng_sbrut E_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	lincom 1/9*(E_2011+E_2012+E_2013+E_2014+E_2015+E_2016+E_2017+E_2018+E_2019)
	local c = r(estimate)
	local ub = r(ub)
	local lb = r(lb)
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "E"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+1
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		gen ub99 = coef + invttail(N, 0.005)*stder
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen ub90 = coef + invttail(N, 0.05)*stder
		gen lb99 = coef - invttail(N, 0.005)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		gen lb90 = coef - invttail(N, 0.05)*stder
		gen lb = `lb'
		gen ub = `ub'
		gen c = `c'
		sort year
		twoway (rcap ub95 lb95 year , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy)) (rarea ub lb year if year >= 2011, color(red%15)) (line c year if year >= 2011, lcol(red%50) lpat(dash)), legend(off) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB2a.pdf", as(pdf) replace
		
restore

	// Figure B2b	

	reghdfe ln_eng_sbrut E_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year ape_n#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	lincom 1/9*(E_2011+E_2012+E_2013+E_2014+E_2015+E_2016+E_2017+E_2018+E_2019)
	local c = r(estimate)
	local ub = r(ub)
	local lb = r(lb)
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "E"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+1
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		gen ub99 = coef + invttail(N, 0.005)*stder
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen ub90 = coef + invttail(N, 0.05)*stder
		gen lb99 = coef - invttail(N, 0.005)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		gen lb90 = coef - invttail(N, 0.05)*stder
		gen lb = `lb'
		gen ub = `ub'
		gen c = `c'
		sort year
		twoway (rcap ub95 lb95 year , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy)) (rarea ub lb year if year >= 2011, color(red%15)) (line c year if year >= 2011, lcol(red%50) lpat(dash)), legend(off) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB2b.pdf", as(pdf) replace
restore		
	
	
	// Figure B2c	
		reghdfe ln_eng_hrs E_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year)  cl(mlabexid)	
		
		
		regsave using "${tmp}/tmp", replace
		lincom 1/9*(E_2011+E_2012+E_2013+E_2014+E_2015+E_2016+E_2017+E_2018+E_2019)
		local c = r(estimate)
		local ub = r(ub)
		local lb = r(lb)
		
	preserve
		use "${tmp}/tmp", clear
			keep if substr(var, 1, 1) == "E"
			split var, p(_)
			destring var2, replace force
			rename var2 year
			local N = _N+1
			set obs `N'
			replace year = ${NORMYEAR} if year == .
			replace coef = 0 if coef == .
			replace stderr = 0 if stderr == .
			replace N = N[1] if N == .
			gen ub99 = coef + invttail(N, 0.005)*stder
			gen ub95 = coef + invttail(N, 0.025)*stder
			gen ub90 = coef + invttail(N, 0.05)*stder
			gen lb99 = coef - invttail(N, 0.005)*stder
			gen lb95 = coef - invttail(N, 0.025)*stder
			gen lb90 = coef - invttail(N, 0.05)*stder
			gen lb = `lb'
			gen ub = `ub'
			gen c = `c'
			sort year
			twoway (rcap ub95 lb95 year , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy)) (rarea ub lb year if year >= 2011, color(red%15)) (line c year if year >= 2011, lcol(red%50) lpat(dash)), legend(off) yscale(fextend) xscale(fextend) ///
			xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
			graph export  "${outpath}/FigureB2c.pdf", as(pdf) replace
			
	restore

	// Figure B2d	

	reghdfe ln_eng_hrs E_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year ape_n#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	lincom 1/9*(E_2011+E_2012+E_2013+E_2014+E_2015+E_2016+E_2017+E_2018+E_2019)
	local c = r(estimate)
	local ub = r(ub)
	local lb = r(lb)
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "E"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+1
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		gen ub99 = coef + invttail(N, 0.005)*stder
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen ub90 = coef + invttail(N, 0.05)*stder
		gen lb99 = coef - invttail(N, 0.005)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		gen lb90 = coef - invttail(N, 0.05)*stder
		gen lb = `lb'
		gen ub = `ub'
		gen c = `c'
		sort year
		twoway (rcap ub95 lb95 year , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy)) (rarea ub lb year if year >= 2011, color(red%15)) (line c year if year >= 2011, lcol(red%50) lpat(dash)), legend(off) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB2d.pdf", as(pdf) replace
restore		
	
	
	// Figure B2e	
		reghdfe ln_eng_hrwg E_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year)  cl(mlabexid)	
		
		
		regsave using "${tmp}/tmp", replace
		lincom 1/9*(E_2011+E_2012+E_2013+E_2014+E_2015+E_2016+E_2017+E_2018+E_2019)
		local c = r(estimate)
		local ub = r(ub)
		local lb = r(lb)
		
	preserve
		use "${tmp}/tmp", clear
			keep if substr(var, 1, 1) == "E"
			split var, p(_)
			destring var2, replace force
			rename var2 year
			local N = _N+1
			set obs `N'
			replace year = ${NORMYEAR} if year == .
			replace coef = 0 if coef == .
			replace stderr = 0 if stderr == .
			replace N = N[1] if N == .
			gen ub99 = coef + invttail(N, 0.005)*stder
			gen ub95 = coef + invttail(N, 0.025)*stder
			gen ub90 = coef + invttail(N, 0.05)*stder
			gen lb99 = coef - invttail(N, 0.005)*stder
			gen lb95 = coef - invttail(N, 0.025)*stder
			gen lb90 = coef - invttail(N, 0.05)*stder
			gen lb = `lb'
			gen ub = `ub'
			gen c = `c'
			sort year
			twoway (rcap ub95 lb95 year , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy)) (rarea ub lb year if year >= 2011, color(red%15)) (line c year if year >= 2011, lcol(red%50) lpat(dash)), legend(off) yscale(fextend) xscale(fextend) ///
			xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
			graph export  "${outpath}/FigureB2e.pdf", as(pdf) replace
			
	restore

	// Figure B2f	

	reghdfe ln_eng_hrwg E_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year ape_n#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	lincom 1/9*(E_2011+E_2012+E_2013+E_2014+E_2015+E_2016+E_2017+E_2018+E_2019)
	local c = r(estimate)
	local ub = r(ub)
	local lb = r(lb)
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "E"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+1
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		gen ub99 = coef + invttail(N, 0.005)*stder
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen ub90 = coef + invttail(N, 0.05)*stder
		gen lb99 = coef - invttail(N, 0.005)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		gen lb90 = coef - invttail(N, 0.05)*stder
		gen lb = `lb'
		gen ub = `ub'
		gen c = `c'
		sort year
		twoway (rcap ub95 lb95 year , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy)) (rarea ub lb year if year >= 2011, color(red%15)) (line c year if year >= 2011, lcol(red%50) lpat(dash)), legend(off) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB2f.pdf", as(pdf) replace
restore		
			